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Abstract. This report concerns the inverse problem of estimating a spatially 
dependent coefficient of a partial differential equation from observations of the 
solution at the boundary. Such a problem can be formulated as an optimal 
control problem with the coefficient as the control variable and the solution 
as state variable. The heat or the wave equation is here considered as state 
equation. It is well known that such inverse problems are ill-posed and need 
to be regularized. The powerful Hamilton-Jacobi theory is used to construct 
a simple and general method where the first step is to analytically regularize 
the Hamiltonian; next its Hamiltonian system, a system of nonlinear partial 
differential equations, is solved with the Newton method and a sparse Jacobian. 
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1. Introduction 

In this paper we study the inverse problem to determine a spatially dependent 
coefficient a of a partial differential equation from partial knowledge of the forward 
solution u. In particular, we seek the diffusion coefficient in the heat equation 
and the wave speed coefficient in the wave equation. Inverse problems arise in 
many applications such as inverse scattering, impedance tomography and topology 
optimization, see e.g. [TJ |3j [fH [14] , and share the property that they are ill posed 
i.e. given data u there may not exist a corresponding coefficient a, and if it exists 
it may not be unique nor depend continuously on u. To be able to determine a 
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the problem thus needs to be regularized such that it becomes well posed. The 
method used here to regularize and to solve the inverse problem is based on the 
work [7J |5J 1151 116) where the inverse problem is formulated as an optimal control 
problem and the corresponding Hamilton-Jacobi equation is used to construct a 
regularization, to obtain convergence results, and to finally solve the regularized 
problem by using the method of characteristics i.e. to solve the corresponding 
Hamiltonian system. 

The paper is stuctured as follows: In Section [2] the general theory of optimal 
control of partial differential equations and Hamilton- Jacobi-Bcllman is presented. 
In Section [3] the idea of how to optimally control the heat equation is discussed 
together with numerical examples, and in Section[4]thc control of the wave equation 
is treated. 



2. Optimal Control and Dynamic Programming 

Consider a differential equation constrained minimization problem with solution 
f : ft x [0,T] -> V, f = <p(x,t) and control a : ft x [0, T] -> B := W, a = a(x,t) 
for an open domain ft, some Hilbert space V on ft, and closed bounded set Bel: 

min / h(f, a) dt + g(f T ), 

(1) a:nx[0,T]^Bj 

ft = f 

with ip T :— </?(•, T) and given initial value ip° = f(-,0). Here, <p t denotes the partial 
derivative with respect to time, / : V x W — » V is the flux, and h : V x W — > R, 
g : V M. are given functions. 

This optimal control problem can be solved either directly using constrained 
minimization or by dynamic programming. The Lagrangian becomes 

L(ip, A, a) := / (A, /(</?, a) — <p t ) + h(<p, a) dt, 



ii 



with Lagrange multiplier A : ft x [0, T] — > V, and the constrained minimization 
method is based on the Pontryagin method 

ft = f(f,cr), 

(2) *t = -(\U{v,<*)) + K( ( P> a )> 

<j(-, t) € argmin{(A, f(<p, a)) + h(ip, a)}. 

with given initial value ip°, final value A T :— A(-,T) = g v (f T ), and where 
h v denotes the Gateaux derivatives with respect to f and (v,w) is the duality 
pairing on V, which reduces to the L 2 (ft) inner product if v,w G L 2 (ft). For a 
differentiable Lagrangian that is convex in a the Pontryagin principle coincides 
with the Lagrangian formulation for a constrained interior minimum 

ft = f(f,<j) 

A* = -(X,f v (f,a)) + h (p (f,a) 
( ' = (X,f <T (f,a))+h (T (f,a), 

aeB, 

but in general ^ and ^ may have different solutions <p, A, a although both describe 
necessary conditions for a minimizer to (|T|). If an explicit minimizcr in (|2]) can 
be found the Pontryagin principle gives additional information about the control. 
Pontryagin's minimum principle can also be written as a Hamiltonian system, see 



SYMPLECTIC RECONSTRUCTION OF DATA FOR HEAT AND WAVE EQUATIONS 3 



0, 

ft = H x (ip,X) 
Xt = -H v (<p,X) 

with ip° given, A T = g v (<p T ), and the Hamiltonian H : V x V — > K defined as 

(5) H(X,(p) := min {(X,f(tp,a)) + h(<p,a)}. 

The alternative dynamic programming method is based on the value function 
U : V x [0, T] -> M, 

CT(^,t):= inf {/ /i(^a) di + 5 (^ T ) ^ = /(p, a), ^(-, r) = € V 

which solves the nonlinear Hamilton-Jacobi-Bellman equation 

(6) d t U(<p,t) + H (U+ fat) ,<f>) =0, U(4>, T) = 5(0), 

with Hamiltonian defined as in ([5]). Note that solving the Hamiltonian system 
Q is the method of characteristics for the Hamilton- Jacobi equation with 
X(x,t) = U,p(ip(x,t),t). In general, the value function is however not everywhere 
diffcrcntiable and the multiplier A becomes ill defined in a classical sense. 

The Hamilton- Jacobi formulation ^ has the advantages that there is a com- 
plete well-posedness theory for Hamilton- Jacobi equations, based on non-differential 
viscosity solutions, see [5], and it finds a global minimum. However, (JsJ) is not com- 
putationally feasible for problems in high dimension, such as the case where ip is 
an approximation of a solution to a partial differential equation. The Hamiltonian 
form Q has the advantage that it is computationally feasible but the drawbacks 
are that it only focuses on local minima and that the Hamiltonian (|5| in general 
only is Lipschitz continuous, even if /, g and /; are smooth, which means that the 
optimal control depends discontinuously on (A, ip) and Q becomes undefined where 
the Hamiltonian is not differentiable. 

In the following sections we will use a regularized version of Q to iteratively 
solve the nonlinear constrained optimization problem ([lj. 

3. Parameter Reconstruction for the Heat Equation 

A distributed parameter reconstuction problem for the heat equation is to find 
a heat conductivity (the control) e.g. a : Cl x [0, T] — ► [er_,er + ], a = a(x,t), 
< <j_ < er+, and a temperature distribution (the state) u : Ct x [0, T] — > V, 
u = u{x,t) that satifies the heat equation 

u t = div(crVu), in x (0,T], 

(7) crVu ■ n = j, on dn x (0, T], 

u = 0, on O x {t = 0}, 

such that the error functional 

(8) / / (u~u*) 2 dsdt, 

Jo Jan 

is minimized. The function u* = u*(x,i) often represents physical measurements 
contaminated by some noise, e.g. u*(x,t) = u trU e{x,t)l + w(x,t) where w is a 
noise term and u trU e satisfies the above heat equation for some unknown parameter 
a true, and in practice the control is only spatially dependent, a trU e = &true{x). The 
primary goal is thus to determine the unknown diffusion coefficient a true and the 
method to do so is to minimize the objective functional ([8]). 

Inverse problems like (|7| , ([8| are in general ill-posed due to one or more of the 
following reasons: 
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(1) There exists no minimizer (it, er), something that may occur with noisy 
data. Given unperturbed data u* corresponding to a true , it is evident that 
there exists a minimizer to ^ , 

(2) The minimizer is not unique, e.g. although it may be possible to find an 
optimal state that minimizes ([8]), u and a may not be unique in ft. 

(3) The solution (u, a), and particularly the control a, depends discontinuously 
on data u* . 

A simple and common way to impose well-posedness to many inverse problems is to 
add a Tikhonov regularization of the form £||o"|||,a(j2x(o T)\ ^ or e > 5 to the objective 
functional ^ , see [TUJ HZJ [TJ [T3] . Using the Pontryagin principle presented in the 
previous section we will in Section 3.2 regularize the inverse problem Q, ^ in a 



way that is comparable to a Tikhonov regularization. 

Formulated as an optimal control problem the most natural assumption on the 
control a is that it is dependent on both time and space but as we will see in Section 



3.3 it is also possible to let a — <r(x), a = a(t), or even let a be constant in time 



and space. 

3.1. The Hamiltonian System. Following Section [2] the Hamiltonian associated 
to the optimal control problem ([7]) and ^ is 

H(u,q,t) := min / (u — u*) 2 ds+ / div(erVit)g dx 

cr:fi->[o-_,CT + ] Jqq Jq 



(9) 



/ (u — u*) 2 + jq ds + min / — uV« • Vg da; 

JdQ. CT:n->[CT_,<7+] Jq 

/ (u — u*) 2 + jq ds — / max {ctVu • Vq} dx. 

JdQ Jf2 ^e[<T_,<r + ] 



ion 

s v ' 

t)(Vu-Vq) 

and the Hamiltonian system, in strong form, then becomes 

u t = div(aVu), inf2x(0,T], 
dVu-n = j, on9Ox(0,T], 

u = 0, on Q x ft = 01, 

(10) ' , 

-q t = div(aVg), in fl x (0,T], 

aVq-n = 2(u-w*), on ffi x (0,T], 
q = 0, on ft x {i = T}, 

with 

(11) 5- := f)'(Vu- Vg). 

It is here evident that the Hamiltonian only is Lipschitz continuous and the control 
a is a bang-bang type control and depends discontinuously on the solutions (u,q), 
see Figure [T| From the optimality conditions ^ an optimal solution has to satisfy 



Vit ■ Vg = and (10) is thus undefined since f)'(0) is set valued, which calls for a 
regularization. 



3.2. Regularization. A simple regularization of the Hamiltonian system ( 10 1, and 
consequently of the Hamiltonian (J9j, is to approximate f)' with the parabolic func- 
tion 

(12) ft(Vu • Vg) := + tanh(^V U • Vg), 
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for some small 5 > 0, see Figure [T] This regularization can be compared with a 
classic Tikhonov regularization where a small L 2 -penalty of the control is added to 
the objective function (tSJ) , i.e. to minimize 



(13) 



(u~u*) 2 ds dt + 5 { [a 2 dxdt. 
o Jan Jo Jn 



Minimizing (13) under the constraint ^ will lead to a C 2 -Hamiltonian with 

H(u,q,t) — / (it — u*) 2 + jq ds — / max { cr( Vu ■ Vq — <5er)} dx, 
Jan Jn o-ek- ,<r+] 



which can be seen in Figure [T] 



tlTihhonoviVu-Vq) 




2 
1.8 
1.6 
1.4 
1.2 



Figure 1. The functions f) (solid line), t)s (dashed line), fyrikhonov 
(dash-dotted line) to the left and their derivatives to the right. 



Another way to describe the simple regularization ( 12 1 is to see what kind of 



penalty on the objective function it corresponds to. We note that the regularized 
Hamiltionian system can be written as 

[ ( [ -u t v - i)' s (\7u ■ X7q)X7u ■ V« dx + [ jv ds] dt = 0, Vw S V, 
Jo \Jn Jan J 

[ ( [ qtV-t)' S (Vu-X7q)X7q-X7v dx + [ 2(u - u*)v ds] dt = 0, Vw S V, 
Jo \Jn Jan J 

or by a redefinition of a 



-Utv — erVu ■ X7v dx + jv ds I dt = 0, Vu S V, 
o \Jn Jan J 

(14) f ( f q t v ~ <rX7q-X7v dx+ [ 2(u - u*)v ds] dt = 0, Vw £ V, 

where cr : [0, T] x £1 — > PF for some Hilbert space W. Let be the primitive function 
of the inverse function of t)' 5 i.e. 
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then it is evident that ( 14 1 can be seen as the first order optimality conditions for 



the problem to minimize 

i-T 



{u-u*y ds dt + 



o Jan 



[ [ #(<r) dx dt, 
Jo Jn 



under the constraint (J7|. In Figure[2]the function S){<j) is compared with a Tikhonov 
regularization of the form 8 (a — of . 



0.35 




FIGURE 2. The function Sj(a) (solid line) compared to the L 2 
penalty function 5(a — a) 2 (dashed line) for 6 = 1, er_ = 1 and 
(7+ = 2. 

It is often beneficial to prevent spacial oscillations of the coefficient by adding 
a penalty on the L 2 -norm of the gradient of the coefficient, i.e. e|| Vcr||^ 2 ( nx ( jy\i 
for e > 0, to the objective function (|8]). For such a penalty the minimization in the 
corresponding Hamiltonian 

(15) H(u,q,t):= min / (u - u*) 2 ds + / div(crVu)g + e\ Vct| 2 dx, 

o- : n-^[cr_,er + ] Jqq Jn 

can not be done explicitly, and instead taking the first variation in a would give 
the system 

u t — div(<7 Vw) , 
— qt = div (crVq) , 
2eAa = -Vu ■ Vq, 
a € [o-_,cr+]. 

which corresponds to the usual first order optimality conditions for the Lagrangian. 
How to treat different penalties on the control in an optimal control setting is 
discussed in Section 13^41 



3.3. Time Independent Control. To study the case when the control a is in- 
dependent of time we first assume that it not only is independent of time but also 
depends on an auxilliary variable z, i.e. a : x [0,T] — > [er_,cr + ], a = a(x,z). 
For a moment we also assume that u : Cl x [0, T] x [0, T] — * V, u = u(x, t, z), but 
with the same measurements as in ([8]) . If we treat z as the time and t as a spacial 
variable we can define the optimal control problem 

i-f e T 



(16) 



mm 

:fix [0,f] 



in — I j I (u — u*) 2 ds dt dz, 

^[<r_,o- + ] T Jo Jo Jan 
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where the state u satisfies the partial differential equation 
1 



— ( div(erVii) — u t 



(17) 



crVii ■ n = j, 
u = Q, 

U = Uq, 



in fl x (0,T) x (0,T], 

on dfl x (0,T) x (0,T], 
on x {t = 0} x (0,T], 
on n x (0,T) x {z = 0}. 



for some arbitrary initial condition u(x,t,0) — uq. 
The Hamiltonian for (16 1, (17) is 

H(u,q,z):= min — / j (u—u*) 

^div(crVu) — q 



o ./n 



2 ds dt 



dx dt 



(18) 



i / / (u- u*) 2 + jq ds dt - i 



o Jo 



dx dt 











[ max J 


H 


f 







f,(f^Vu-Vq dt) 

and the Hamiltonian system is given by 

^(div(F)'Vu) -inj, 



(19) 



- j> 
t)'V« • n = j, 

u = 0, 

U = MO 
1 



-<Z, = ,~(div(f)'V3) + g t 

fj'V<z-n = 2(it-u*), 
9 = 0, 
9 = 0, 



in ft x (0,T) x (0,T], 

on dfl x (0,T) x (0,f], 
on 0, x {£ = 0} x (0,f], 
on ft x (0,T) x {z = 0}, 

in ft x (0,T) x (0,f], 

on dft x (0,T) x (0,T], 
on ft x {i = T} x (0,T], 
on ft x (0,T) x {z = f}. 



Under the assumption that the solutions u and q in ( 19 I are asymptotically sta- 
tionary as T — > oo, the Hamiltonian system for the problem (17]), (raj, with a time- 
independent control, is given by (10 1 and 



(20) 



Vit ■ Vq dt 



Similarly, the case of a space independent coefficient a — a(t) will lead to 



and for the case where a is constant 



Vm • Vq dx I , 



Vw • Vg da; dt 
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3.4. Penalty on the Control. If wc want to reconstruct a time independent 
control it can be beneficial to put a penalty on at, i.e. we want to minimize the 
objective functional 



(21) 



F(u,a t ):= [ [ (u- u*) 2 ds dt + e [ [ a\ dx dt, 
Jo Jan Jo Jn 



under the usual constraint ([7]). To do this the optimal control problem has to be 
reformulated such that a is a state variable and the control is defined as z := tx t (x, t), 
z : SI x [0, T] — » The optimal control problem is thus to find a control z 

and state variables u and a such that F(u, z) is minimized and the system 

u t = div(crV-u), in SI x (0,T], 

a t = z in fix (0,7*], 

crVu-n = j, on 3fix(0,T], 

u = 0, on x {t = 0}, 

cr = cr > 0, on Cl x {t = 0}. 
is satisfied. The Hamiltonian becomes 

H(u,q,a,X,t) := min / (tt — u*) 2 ds+ / div(crVw)(7 + zA + ez 2 dx 

z:n-.[ 2 _ : z + ] Jqq Jf2 

= / (u — u*) 2 + jq ds — / aX7u-Vqdx 
Jan Jn 

min {z(ez + A)} dx, 



and the corresponding Hamiltonian system is 



Ut 


= div(crVu) , 


in fl x (0,T], 


<rt 


= f)'(A) 


in Q x (0,T], 


(tVm • n 


= ij 


on 90 x (Q, T] , 


u 




on Q x {i = 0}, 


a 


= <Jo > 0, 


on x {i = 0}. 


-Qt 


= div(crV<7) , 


in SI x (Q, T] , 


-A t 


= -Vit • Vq 


in ft x (0,T], 


<tV<7 • n 


= 2(u-u*), 


on dfl x (0,T], 


q 


= 0, 


on SI x {t = T}, 


A 


= 0, 


on SI x {t = T} : 



which is equivalent to (10) with 



lo \ly 



o-:=<t +/ t)' / -(Vm- Wq){x,z) dz } dy 



Note, since we no longer have a constraint a > 0, the bound z_ has to be carefully 
chosen to ensure well-posedness of the forward problem. 

In a similar fashion as for penalizing temporal variations of the control it is 
also possible to penalize spacial variations, as was briefly mentioned in Section [3. 2| 
where the objective was to minimize F(u, |Vct|) under the constraint which 



leads to the Hamiltonian (151. To be able to explicitly find the minimum in the 
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Hamiltonian we once again let tr act as a state variable, introduce the control z and 
the dynamics 

, s ot= Z ~ |W|2 , infix (0,T], 

(22) 7 V ' J ' 

cr = (T > 0, in fi x {t = 0}, 

for 7 > 0. The slightly perturbed control problem is now to minimize the objective 
function F(u,z) such that ([7| and (22 1 holds, which leads to the Hamiltonian 

f f *^ I ^ ^ I ^ 

H(u, q, a, A, t) :— min / (u— u*) 2 ds + / div (crVw) q + A \- ez dx 

z-.n^[z-,z + ] J m J n ' 7 

= / (u - u*) 2 + jq ds - [ aVu ■ Vq + dx 
Jan Jn 7 

+ [ min {z(e H — )} da;, 

J n z:n^[z-,z + ] 7 



and the Hamiltonian system 



f)(A) 

u t = div(crVw) , 

= m - , 

7 

—q t = div(crVg), 

Act 

-A t = Vu- Vg-2A — . 

7 

3.5. Numerical Approximation and Symplectic Methods. Let V C V := 

i? 1 (fi) be the finite clement subspace of piecewise linear functions defined on a 
triangulation of fi, which implies that our optimal control problems in the previ- 
ous sections are approximated by optimal control problems for ordinary differential 
equations. We also let the functions i)s,H s and hg denote the regularized counter- 



parts to i),H and h. The regularized version of f) is given by (12 1 from which the 



definition of H follows. The regularized function hg can be derived from H by 
hg := H 5 — (A, H s x ) and a regularized version of / can be defined as fg := H 5 X . 

Now, introduce the uniform partition {ti — ki}fL , k — T/N of the time interval 
[0, T], and the corresponding finite clement approximations at each time step ip n := 
if(t n ),X n :— X(t n ). Also define a discrete regularized version U : V x [0, T] — > M of 
the value function ([2|, 

, N-l 

U(4>,t m ) := min \g{(p N ) + k S~] hg((p n ,X n +i) 



n—m 



where tp n and A„ satisfy a symplectic scheme, e.g. the symplectic forward Eulcr 
method 

ip n+ i - tp n = kH s x (ip n , A„ + i), for n = m, . . . , N - 1 given ip m = (f>, 

A„ - A n+ i = kH Jip n , A„+i), for n = m, . . . , N - 1 given X N = g v (<fiN)- 



Symplecticity here means that Uu,(ip n ,t n ) = A„, i.e. the gradient of the discrete 
value function coincides with the discrete dual A„, and given that \H — H s \ = 0(6) 
it can be shown that for symplectic one-step schemes 

JV-l 

0(k), 



U(<po, to) - g(ipN) ~ky~] hg((p n , A„+i) 
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for 5 ~ k, see [T5]. It is thus essential to use a symplectic time discretization of the 
regularized Hamiltonian system 

in order to have convergence in the value function. 

Some examples of other symplectic schemes are the the backward Euler method 

(pn+i - Vn = kH s x ((p n+1 ,X n ), for n = 0, . . . , N - 1 given ip , 
(24) s 

A„ - A n+ i = kH v (if n+1 ,X n ), for n = 0, . . . , N - 1 given X N , 

and the implicit midpoint method 
(25) 

, ttS ( Pn + ¥>n+l X n + X n+ i \ at t ■ 

Vn+i ~ <p n = kH x , , for n = 0, . . . , N - 1 given ip , 



\ \ i ttS ( V>n + X n + X n+ i . . 

X n - X n+1 = kH^ , , for n = 0, . . . , N - 1 given X N . 



See |12j for a thorough description of symplectic methods. 

3.6. The Newton Method. To solve the coupled nonlinear symplectic schemes 
( 23 1-( 25 ) above, it is tempting to propose fix-point schemes that partly removes the 
coupling between the forward and bacward equation, e.g. by iterating separately in 
ip and A. Such methods has the advantage that existing partial differential equation 
solvers can be used to efficiently solve the forward and backward problems in each 
iteration, but the disadvantage is that the convergence to an optimal solution tends 
to be slow, and also dependent on the discretisation. A more suitable strategy is 
to use information of the Hessian of H 5 ; e.g. Quasi-Newton methods, or since the 
Hessian in our case can be found explicitly and is sparse, the Newton method itself. 

For the Hamiltonian system (10 1 with a := given by (12 1 the symplectic 
backward Euler can be written as 

F n (w) = 0, G n (w) = 0, n = Q,...,iV- 1, Vw e V 



where 



F n (w) := / (u n +i - u n )w + fcf^(Vu n+ i • Vq n )Vu n+ i ■ Vw dx 
Jn 

- / kj n+1 w ds, 
Jan 

G n (w) := / (q n - q„+i)w + if)J(Vfi n+ i • Vq n )Vq n ■ Vw dx 
Jn 

2k(u n+1 - u* n+1 )w ds, 



(26) 



I an 

and wo = qN — 0. Given an initial guess u[0], q[0] the (damped) Newton method 
yields that 

u[i + 1] = u[i] — ail, 
q[i + 1] = q[i] - aq, 

where a G (0, 1] and, for each iteration, the updates u and q solve a linear system 
of the form 



(27) 



K n K u \ ( u \ ( f 
K 21 )\q \g 
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where 

u = ( ui ... u N ) , q = ( q ... giv-i ) , 

/ = ( Fq ... F/v-i ) , g = ( G ... Gjv-i ) ■ 

The matrix Kn is a bi-diagonal block matrix with M + Si for i — 0, . . . ,N — Ion 
the diagonal and — M on the sub-diagonal, where M denotes the mass matrix 



ww dx, 



and 



S n ■= / Ws(Vun+i ■ V<j„)V<j n • Vro Vu„+i • Vw dx 
Jq 



+ / kt)' s (Vu n +i ■ Vq n )Vw ■ Ww dx. 
Jn 

for w,w £ V. The matrices K 12, i^i are symmetric block-diagonal matrices with 

/ kt)g(Vu n +i ■ V9„)Vu„ + i • Ww Vw„ + i • Vw da;, 
Jq 

and 



for n = 0. 



/ fcf)j(V«„+i • Vg„)Vg„ • X7w Vq„ • Vw dx — 2kww ds, 
Jn Jon 

1 on the the diagonals, respectively. 



If we repartition the block 2x2 linear system (27) to 
(28) 



K21 



we see that it is a generalized saddle point system [3] with symmetric matrices 
K2i,Ki2, and Kf x ^ 0, K21 ^ 0. However, unlike saddle point problems aris- 
ing from e.g. the steady-state Navier-Stokes equations or from the Karush-Kuhn- 
Tucker optimality conditions for equality constrained minimization problems, both 
K12 and K21 may here be indefinite and singular. 



Since ( 27 ) and ( 28 1 are increasingly ill-conditioned with respect to reduction in 



mesh size, step size and regularization, the success of iterative algorithms like Krylov 
sub-space methods will depend heavily on the choice of preconditioner. Standard 
algebraic preconditioners like incomplete LU-factorization are often unsuitable for 
saddle-point problems due to the indefiniteness and lack of diagonal dominance, so 
the preconditioner must be tailored for the specific problem at hand. One popular 
approach for PDE-constrained optimization problems is to base the preconditioner 
on the solution from a reduced approximated problem where the Schur complement 
is replaced by an approximation e.g. by quasi-newton methods, see [5]. 



In our case we use the GMRES method to solve the non-symmetric system ( 27 ) 



and base our preconditioner on the approximate solution of a simple blockwise 
Gauss-Seidel method i.e. to start with a guess q° and iteratively solve 

' l+1 = f-K l2 q\ 



(29) 



KM 



g - K 2 \u 



i+l 



which works well for large regularizations i.e. when t)g is small and the diagonal 



blocks of (27 1 are dominant. Also, each iteration with this method only requires 



one forward and one backward solve in time of a modified heat equation so the com- 
putational work for one iteration is concentrated to solving N — 1 smaller systems 
with system matrices (M + Si). In practice, the Gauss-Seidel method will break 
down for small regularizations but for our problems (and discretizations) only one 
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iteration with ( 29 1 turns out to be a fairly good approximation to use as precondi- 



tioner. Note that for q° = 0, one Gauss-Seidel iteration is the same as solving (27) 
with the approximation K\ 2 = 0. 

Another more elaborate idea is to use a preconditioner based on the solution of 
an approximated Schur complement system 

K n K l2 \(u\ = ( g 
S )\q ) \f- K 12 K^g 

where S is an approximation of the Schur complement 

K[ x - K 21 K^K 12 . 

which essentially is to find a good approximation of the lower triangular block 
matrix . 

Although solution algorithms for saddle point systems on the symmetric form 



(28 I are extensively treated in the litterature, see jj] for an overview, we here favour 



the non-symmetric form ( 27 1 , since a Schur complement reduction of ( 28 ) means 



to find an approximation to the Schur complement 

which since K 2 \ here can be singular, is unavailable. One way around this obstacle 



is to rewrite (28 1 by e.g. the augmented Lagrangian method which leads to a 
symmetric invertible Schur complement but where the physical meaning of the 
original system, on PDE level, is partially lost. 



If a direct solver is used for the Newton system it is appropriate to reorder ( 27 ) 
such that the solution vector and right hand side contains time steps in increasing 
order, which leads to a banded Jacobian with band-width of the same order as the 
number of spacial degrees of freedom. 

Our computations were implemented MATLAB (for the one dimensional exam- 
ples), and in DOLFIN [13], the C++/Python interface of the finite element solver 
environment FEniCS [TT] (for the two dimensional examples) . Piecewise linear ba- 
sis functions were used for the finite element subspace V, and in all examples the 
solution u, q was first calculated for a large regularization which was succcsively 
reduced such that the solution from the previous regularization served as starting 
guess for a smaller regularization. 



For the two dimensional examples the sadde-point system ( 27 I was solved with 
the PETSc implementation of GMRES (used by DOLFIN) with preconditioning 
based on the solution from one iteration of blockwise Gauss-Seidel method. For 
the one dimensional examples a direct solver was used. The number of iterations 
for GMRES with the Gauss-Seidel preconditioner seems to be relatively insensitive 
with respect to temporal and spacial discretization but still highly sensitive to the 
regularization in our examples. 

To give a time independent approximation <j(x) of the time dependent control 
a{x,t), approximated by a := r/ 5 (Vit- Vg) where u,q are solutions to the Hamilton- 



ian system (10 1, three different types of averaging were tested as post-processing: 



(1) Let the time independent control be defined by the Hamiltonian (18 1, i.e 



(30) °"' =t, ' S {J \7u-X7qdt 

(2) Let the time independent control be the average of the time dependent 
control, i.e. 

(31) a:= hj ft(V«-V?)dt. 
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(3) Let the time independent control be the weighted average 

fm , / T ^(V M -V g )|V M -V g | dt 
(62) a := = , 

filVu-Vq] dt 

of the time dependent control f)J>(Vu • Vg). 

The weighted average turned out to be the most successful aproximation and can 
be explained by first extending the Hamiltonian ^ to also depend on the artifical 
variable z as in Section [3.31 

H(u,q,z):=— j j (u — u*) 2 + jq ds dt j j u t q dx dt 

T Jo Jan T Jo Jn 

- i [ [ f)'(Vit • Vg)Vu ■ Vg dt dx, 
T Jn Jo 

where f)'(Vu ■ Vg)Vu ■ Vg = f)(Vit • Vg) by definition. For the problem with a time 
independent control we now seek an approximation of the Hamiltonian ( 18 1 of the 
form 

H(u,q,z):— — j f (u — u*) 2 + jq ds dt ~ f f u t q dx dt 

T Jo Jan T Jo Jn 

- i / /(Vu • V<j) / Vu • Vg dt dx, 
T Jn Jo 

that best approximates H, i.e. 

fm v , Jo WJu ■ Vg)Vu ■ Vg dt 

f{Vu ■ Vg) := -jt— — — . 

J Q Vu • Vg dt 

In Figure [3j one dimensional reconstructions from three sets of simulated data 
u* , generated from a time independent conductivity <Jt rue , are compared: 

(1) Data calculated with the same discretization as u and q. 

(2) Different discretisations used for data and solutions. 

(3) Different discretisations used for data and solutions and with noise in the 
data u*. 

The last set is the most realistic one since for true experimental data of u* it is 
inevitable to not only have measurement noise but also a systematic error from 
the numerical method. To simulate noise the discrete solution u* was multiplied 
componentwise by independent standard normal distributed stochastic variables 
Eij according to u*(xi,tj)(l + Tjeij), where rj denotes the percentage of noise. It 
is notable that the systematic error from using different meshes can have a much 
bigger effect on the solutions than additional noise, which can be observed from the 
dual solution q in Figure [3] 

In Figure [4] the time independent post-processing of the time dependent recon- 
struction can be found. It is here evident that the weighted average (32 1 performs 



better than ( [31] ) , but since the reconstruction is highly dependent on the given 
boundary condition, see Figure [5] for comparison, there are situations where the 
different post-processing techniques perform equally well. It would of course be 
optimal to use the knowledge that Otrue. is independent of time in the calculations, 
i.e. to use the Newtion system for (10 1 with time independent-control (20). This 
would however lead to a dense Jacobian. 

Note that in the examples the limits <j_ , er + were chosen to be the biggest and 
smallest values of <J true . In our experience the Pontryagin method is not well 
suited for reconstruction of values between cr_ and <r+ if there is noise or other 
measurement errors present in data. 
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Figure [6] shows two-dimensional reconstructions of two different time indepen- 
dent conductivities. Unlike the one-dimensional example the quality of the recon- 
struction here deteriorates quickly as the distance to the measurement locations is 
increased. 




Figure 3. ID reconstruction of a trU e — 0.75 — 0.5tanh(20a; — 10) 
with S = 1CP 6 , <7_ = 0.5, (7+ = 1, measurements on both bound- 
aries and Neumann boundary condition au x (0,t) — ~au x (l,t) — 
sin(4t) for t < 0.5 and elsewhere. The plot shows, from top to 
bottom, it, q, t)' 6 and the objective function \\u — u*||l2(3q x [o i t])- 
In all cases it, q was calculated with 50 steps in space and time. 
In the left column, the data u* was generated by solving the heat 
equation with 50 steps in time and space and conductivity a true, 
while 200 steps in time and space was used in the middle and right 
columns. In the right column 10% noise was also added to u* . 



4. Reconstruction from the Wave Equation 

In this section the goal is to determine the wave speed for a scalar acoustic wave 
equation: Given measured data u* , find the state u : f2 x [0, T) — > V, u — u(x,t) 



SYMPLECTIC RECONSTRUCTION OF DATA FOR HEAT AND WAVE EQUATIONS 15 



D.95 ■ 

0.9 , ' 
0.85 ■ 




°" X 



Figure 4. The time independent post-processed conductivity for 
the ID reconstructions in Figure [3j The true control a true is indi- 
cated by a solid line and the averaged controls (30 1, (31) and (32 1 



are indicated by dotted, dash-dotted and dashed lines, respectively. 





x 





x - 



x 1 



FIGURE 5. ID reconstruction with data as in Figure [3] and [J] but 
with Neumann boundary condition au x (0,t) — —au x (l,t) = 1. 




and a control e.g. a : x [0,T] — » [cr_,o+], u = <r(x,t) where < er_ < er+, that 
solves the partial differential equation 

u u = div(crVu), 

(33) crVu ■ n = j, 

u = Ut = 0, 

such that the error functional 



in n x (0,T], 
on dfl x (0,T], 
on!ix{( = 0}, 



(34) / / (u-u*) 2 ds dt, 

JO JdQ. 

is minimized. The control a is here the square of the wave speed of the medium 
and u is the pressure deviation. 



To use the framework of the previous section we note that ( 33 1 can be written 
as the first order system 

v t = div(crVu), inf2x(0,T], 
u t = v, infix(0,T], 
aVu-n = j, on 90 X (0,7*], 

U = v = 0, on n x {t = 0}. 



(35) 
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Figure 6. 2D reconstruction on the unit square with final time 
T = 1 and Neumann boundary condition cr|^ = 1 on dfl x [0,T]. 
The data u* was simulated by solving the forward equation on a 
quasi-uniform mesh with 13000 triangles and 80 time steps while 
the inverse problem was solved on a uniform mesh with 3200 tri- 
angles and 40 time steps. Measurements from the whole boundary 
were used. Top: True conductivity atrue- Middle: Reconstructed 
condictivity for S w 0.002 using the weighted average (32 1. Bot- 
tom: As in middle but for i5 ~ 0.05 and with 5% noise in the 
measurements. 
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4.1. The Hamiltonian System. As in Section 3.1 we have a Hamiltonian asso- 
ciated with the optimal control problem ( 34 1 and ( 35 1 which is defined by 

(u — u*) 2 ds + / div(crVu)g + vp dx 



H 



mm 

er:Q — ► [cr_ ,er+] 



OP. 



(36) 



(u — u*) 2 + jq ds + 



mm 

:ft — >[(T_ ,(T4 



/ (u — u*) 2 + jq ds + / vp — max {crV 

Jan Jsi <re[<r-,(T + ] 



rrVu ■ + up dx 
Vq} dx, 



and the Hamiltonian system becomes 



(37) 



or equivalently 



(38) 



with 



v t ■ 
u t ■ 
ctVu • n : 
u ■ 
-Vt ■ 

-Qt - 
aVq ■ n : 

V - 

u tt ■ 
(tVu • n : 

u ■ 

qtt 
crVq • n : 

q - 



div(crVu) . 

■■ v, 
■ 3, 

v = 0, 
div(CTVg) , 
P, 

: 2(U-U*), 

9 = 0, 

div(crVu) . 

: 3, 

u t = 0, 
div (crVg) , 

: 2(U-U*), 

Qt = 0, 



fj(Vu-Vg) 

in ft X (0,T], 
in SI x (0,T], 
on dCl x (0,T], 
on x {t = 0}, 
in SI x (0,T], 
in x (0,T], 
on 8Q x (0,T], 
on SI x {t = T 7 }, 

in SI x (0,T], 
on 90 x (0,T], 
on SI x {£ = 0}, 

in SI x (0,7*], 
on <9S! x (0,T], 
on Si x {i = T}. 



(7 := f)'(Vu • Vg). 

4.2. Symplecticity for the Wave Equation. As a natural case the symplectic 
methods discussed in 3.5 with ip = (u, v), A = (p, q), can be used to solve the system 
( 37 1 . It is however also possible to use a time-discretization that is symmetric in 
time i.e. 



(39) 



for & n :— t)'(Vu„ • Vg n ) and n — 1, . . . , N — 1. For a given a, constant in time, this 
scheme is the symplectic backward Euler method for the forward wave equation for 
u, which can be written as the Hamiltonian system (35 1 with Hamiltonian 

1 

2 



u n+ i - 2u n + u n -i 


= fc 2 div(<7„Vit n ), 


in SI, 


d n S7u n ■ n 




on 9S1, 


"0 


= ux = 0, 


in SI, 


Qn+i ~ 2<jv» + q n -i 


= /cdiv(o-„V<7„), 


in SI, 


5-nVq n ■ n 


= 2(u n - u*), 


on 9S1, 


<1N 


= gjv-i = 0, 


in SI, 



,(u,v) :-- 



IctVwI 2 + v 2 dx, 
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and the symplectic forward Euler method for the backward wave equation for q. 



To see that that the scheme (39 1 is symplectic for a n :— r/(Vu n • Vq n ) we note 
that a one-step method (tp n ,\ n ) — > (ip n+ i, A„ +1 ) is symplectic if there exists a 



function H (ip ni A„ +1 ) such that (23 1 holds, or equivalently H((p n+ i, A„) such that 
( 24 1 holds, see Remark 4.8 in [TS] or [T^] for details. It thus follows that the one-step 
method 

v n +i - v n = fcdiv(f)'(Vu„ ■ Vg„)Vu„), 
u n+ i - u n = kv n+ i, 
Pn - Pn+i = fcdiv(f)'(VM„ • Vg„)Vg„), 
q n ~ q n +i = kp n+ i, 
corresponds to the symplectic forward Euler method for the Hamiltonian 

H(u n ,q n ,v n+1 ,p n+1 ) := H(u n ,v n+1 ,p n+1 ,q n ) - 2 / v n+1 p n+1 dx, 
' v v ' Jn 



where H is given by ( 36 1 . Since ( 39 1 only is stable for sufficiently small time-steps 
and still requires to solve a complex saddle point system we will use the symplectic 
midpoint method in our experiments. 



4.3. Numerical Examples. Let a := f)^ where ty s is given by ( 12 1 . The symplectic 



on 



midpoint method for the regularized Hamiltonian system (37) can then be written 
as 

F„ 1 H=0, F„ 2 H = 0, G r \H = 0, Gl(w) = 0, 
for n = 0, . . . , N — 1, and Vw £ V, where 

F n( w ) : = jjy n+ i - v n )w + k\)' s (Vu n+ i ■ Vq n+ i)Vu n+ i ■ Vw dx 
kj n+ iw ds, 

(lln+l - u n - kv n+ i)w dx, 

(q n - Qn+i - kp n+ i)w dx. 

G 2 n (w):= / (p n - p n+1 )w + kt)' s (Vu n+ i ■ Vg n+ i)Vg n+ i • Vw dx 
Jo 

- / 2k{u n+ i -u* H )w ds, 
Jon 2 

and Uq — vq = p^ = qN = 0. The index n + \ implies the average of the values 
at n and n + 1, i.e. u n+ i :— \{u n + u n+ i). Taking the variations with respect to 
u,v,p,q gives the Newton system 



( K n 


K 12 





K u \ 








( h\ 


K 21 


K 22 










V 




h 





Q 


^33 


^34 




p 




91 


V Ku, 





^43 


K i4 J 




K <? ) 




\ 92 ) 



(40) 



with increments 

u = (u 1 ... ii N ) T , v = (vi ... i> N ) T , 
P = (Po ■■■ Pn-1 ) T , Q = (qo ■■■ q~N-i ) 7 
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and right hand side 

/l=( ^0 ■■■ F N-1 T > h = { F ■■■ F N-1 ) T > 

9i = ( Gq ... G 1 N _ 1 ) , g 2 — ( Go ... G 2 N _ X ) . 
with submatrices with the following structure: 
• ifn is lower block bi-diagonal with 



(41) 



on its main diagonal for n = 0, . . . , N — 1 and on its sub-diagonal for 
n= 1,...,N- 1. 

• if 44 is upper block bi-diagonal with (41 1 on its diagonal for n = 0, . . . , N— 1 
and on its super-diagonal for n = 0, . . . , N — 2. 

• if 12 = if 21 = ^ 34 — ^ 43 is lower block bi-diagonal with mass matrices M 
on the main diagonal and — M on the subdiagonal. 

• if 22 = if 33 is lower block bi-diagonal with —^f- on the diagonal and the 
sub-diagonal. 

• if 14 is upper block bi-diagonal with 

\ I ^i'(Vti n+ i ■ Vg I1+ i)Vu n+ i ■ X7w Vm„ + i ■ Vw dx, 

J Q 

on its diagonal for n = 0, . . . , N — 1 and on its super-diagonal for n = 
0,...,N-2. 

• if4i is lower block bi-diagonal with 

I [ ^s(V\ + i • V<?„ + i)Vg„ + i • Vic Vg„ + i • Vw da; - / kww ds, 

on its diagonal for n = 0, . . . ,N — 1 and sub-diagonal for n = 1, . . . , N — 1. 
As in the previous section we will solve the Newton system using GMRES and 
an approximate solution as preconditioner, e.g. from the the 2x2 blockwise Gauss- 
Seidel method 



K n u« 


- 1 - 


- Kl2^ 


-1 


= A- 


ifi49 l , 




- 1 - 


- K 22 v l ~ 


-1 


= /a, 




K33P- 


hi . 


h if 34 <T 


-1 


= 




K 43 f- 


hi . 




-1 


= .92 ~ 


K 41 u l+1 



which can be written as 

(if n - K 12 K^K 21 )u i+1 = fi- K 12 K^f 2 - if 449' 

( 42 ) 1 - , 1 , 

(if 44 - K 43 K^K 34 )q l+1 = g 2 - K 43 K^ 9l - K 41 u 



Note that (42 1 is easily solved since inverting if 22 and if 33 only involves the cal- 
culation of M~ l . In fact, the Schur complements ifn — ifi2if 22 1 if2i and ifi4 — 
K 43 K^ K 34 becomes lower and upper block trianglar matrices, respectively, and 
(42 1 can be solved by one forward substitution in time for u l+1 and one backward 
substitution in time for q l+1 . Of course, to save memory the Schur complement 
system (42) should never be formed explicitly. For large regularizations the Schur 
complements can be seen as approximations of the operator —A + du- As for the 
case with the heat equation starting with q° = 0, one iteration with (42 1 is the 
same as solving (|40| with if 14 = 0. 



In Figure [7J a two dimensional example of reconstruction two different speed 
coefficients is shown. The measured data was here simulated by solving the wave 
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equation for a true with the symplectic backward Eulcr method for ( 35 1 , which can 
be written as the second order scheme 

/ (u n +i — 2u„ + u n -.\)w dx = / jw ds ~ / ctVu„ • Vw dx, Vw G V. 
Jn Jan Jn 

Since the wave equation is a conservation law and is reversible in time it is tempting 

to believe that it would be easier to control than the heat equation but there are 

some computational drawbacks: numerical errors are propagated in time and there 

seems to be many local minima. From the approximation t)' s (Vu ■ Vg) in Figure [8] 

it is evident that the time dependent reconstruction varies a lot over time and is 

not a good approximation of the time independent wave coefficient a true- 
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0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



X X 

FIGURE 7. 2D reconstruction using the weighted average ( |32| , fi- 
nal time T = 1.5 and Neumann boundary condition 2sin(47rf) for 
x G [0.4,0.6], t < 0.5 and elsewhere. The data u* was sim- 
ulated by solving the forward equation on a quasi-uniform mesh 
with 3232 triangles and 328 time steps while the inverse problem 
was solved on a uniform mesh with 1250 triangles and 30 time 
steps. Measurements from the whole boundary were used. Top: 
Reconstruction of a tT ue — 0.5 inside the square [0.2,0.5] x [0.5, 0.8] 
and a true — 1 elsewhere, with no noise in data (left) and 10% noise 
in data (right). Bottom: Reconstruction of Otrue — 0.5 inside the 
square [0.35,0.65] x [0,0.3] and a trU e — 1 elsewhere, with no noise 
in data (left) and 10% noise in data (right). 



22 



JESPER CARLSSON 




Figure 8. Measurements u* (top) and f)^(Vu • Vq) (bottom) for 
timesteps 5, 15 and 25. The data here corresponds to the top left 
plot in Figure [7] and u* is interpolated onto the mesh used for the 
calculation of u and q. 



SYMPLECTIC RECONSTRUCTION OF DATA FOR HEAT AND WAVE EQUATIONS 23 



References 

[1] H. T. Banks and K. Kunisch. Estimation techniques for distributed parameter systems, vol- 
ume 1 of Systems & Control: Foundations & Applications. Birkhauser Boston Inc., Boston, 
MA, 1989. 

[2] Emmanuel Nicholas Barron and Robert Jensen. The Pontryagin maximum principle from 

dynamic programming and viscosity solutions to first-order partial differential equations. 

Trans. Amer. Math. Soc, 298(2):635-641, 1986. 
[3] M. P. Bends0e and O. Sigmund. Topology optimization. Springer- Verlag, Berlin, 2003. Theory, 

methods and applications. 
[4] Michcle Benzi, Gene H. Golub, and Jorg Liesen. Numerical solution of saddle point problems. 

Acta Numer., 14:1-137, 2005. 
[5] George Biros and Omar Ghattas. Parallel Lagrangc-Ncwton-Krylov-Schur methods for PDE- 

constraincd optimization. I. The Krylov-Schur solver. SIAM J. Sci. Comput., 27(2):687-713 

(electronic), 2005. 

[6] Liliana Borcea. Electrical impedance tomography. Inverse Problems, 18(6):R99-R136, 2002. 
[7] Jesper Carlsson. Pontryagin approximations for optimal design of clastic structures, preprint, 
2006. 

[8] Jesper Carlsson, Mattias Sandbcrg, and Anders Szepessy. Symplectic pontryagin approxima- 
tions for optimal design. Preprint. 

[9] M. G. Crandall, L. C. Evans, and P.-L. Lions. Some properties of viscosity solutions of 
Hamilton- Jacobi equations. Trans. Amer. Math. Soc, 282(2) :487-502, 1984. 
[10] Heinz W. Engl, Martin Hanke, and Andreas Neubauer. Regularization of inverse problems, 
volume 375 of Mathematics and its Applications. Kluwer Academic Publishers Group, Dor- 
drecht, 1996. 

[11] FEniCS. FEniCS project. URL: urlhttp//www. fenics.org/. 

[12] Ernst Hairer, Christian Lubich, and Gerhard Wanner. Geometric numerical integration, vol- 
ume 31 of Springer Series in Computational Mathematics. Springer- Verlag, Berlin, second 
edition, 2006. Structure-preserving algorithms for ordinary differential equations. 

[13] J. Hoffman, J. Jansson, A. Logg, and G. N. Wells. DOLFIN. URL: url- 
http//www. fenics.org/dolfin/. 

[14] J.-L. Lions. Optimal control of systems governed by partial differential equations. Translated 
from the French by S. K. Mitter. Die Grundlehren der mathematischen Wissenschaften, Band 
170. Springer- Verlag, New York, 1971. 

[15] M. Sandberg and A. Szepessy. Convergence rates of symplectic Pontryagin approximations 
in optimal control theory. M2AN, 40(1), 2006. 

[16] Mattias Sandberg. Convergence rates for numerical approximations of an optimally controlled 
Ginzburg-Landau equation, preprint, 2006. 

[17] Curtis R. Vogel. Computational methods for inverse problems, volume 23 of Frontiers in 
Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, 
PA, 2002. With a foreword by H. T. Banks. 

CSC, Numerical Analysis, Kungl. Tekniska Hogskolan, 100 44 Stockholm, Sweden; 
E-mail address: jesperc@kth.se 



